
***********************************************************************************
* ML weibull d0/d1 allowing for potential interval censoring and HS heterogeneity *
***********************************************************************************

* single solution, since heterogeneity gap is contrained to be positive
* logistic tranformation for probability term built into the log likelihood

program duration, eclass
 args todo b lnf g
 
 local t_0 "$ML_y1" // duration ellapsed before the observation window
 local t_a "$ML_y2" // start of window where transition takes place
 local t_b "$ML_y3" // end of window where transition takes place

 tempvar x a u o // parameters
 mleval `x' = `b', eq(1)        // weibull location
 mleval `a' = `b', eq(2) scalar // weibull shape
 mleval `u' = `b', eq(3) scalar // log of unobserved heterogeneity for type 1 (higher hazard) => u = exp(ln_u)
 mleval `o' = `b', eq(4) scalar // log odds of being type 1 (higher hazard) => p = exp(`o') / (1 + exp(`o') ) 
 
 tempvar lj // individual contribution to the likelihood according to the censoring case and probability to be of a given type
  
 // [A] flow sampled, uncensored
 qui gen double `lj' = log(             (exp(`x'           ) * `a' * (`t_a' ^ (`a' - 1)) * exp(-exp(`x'           ) * (`t_a' ^ `a')))  ///
                           + exp(`o') * (exp(`x' + exp(`u')) * `a' * (`t_a' ^ (`a' - 1)) * exp(-exp(`x' + exp(`u')) * (`t_a' ^ `a')))) ///
                     - log(1 + exp(`o')) ///
 if `t_0' == 0 & `t_a' == `t_b'

 // [B] flow sampled, right censored
 qui replace    `lj' = log(             (exp(-exp(`x'           ) * (`t_a' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_a' ^ `a')))) - log(1 + exp(`o')) ///
 if `t_0' == 0 & `t_a' <  `t_b' & missing(`t_b')

 // [C] flow sampled, interval censored
 qui replace    `lj' = log(             (exp(-exp(`x'           ) * (`t_a' ^ `a')) - exp(-exp(`x'           ) * (`t_b' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_a' ^ `a')) - exp(-exp(`x' + exp(`u')) * (`t_b' ^ `a')))) ///
                     - log(1 + exp(`o')) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' != 0 & !missing(`t_b')

 // [C2] flow sampled, left censored
 qui replace    `lj' = log(             (1 - exp(-exp(`x'           ) * (`t_b' ^ `a')))  ///
                           + exp(`o') * (1 - exp(-exp(`x' + exp(`u')) * (`t_b' ^ `a')))) - log(1 + exp(`o')) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' == 0 & !missing(`t_b')
        
 // [D] stock sampled, uncensored
 qui replace    `lj' = log(             (exp(`x'           ) * `a' * (`t_a' ^ (`a' - 1)) * exp(-exp(`x'           ) * (`t_a' ^ `a')))  ///
                           + exp(`o') * (exp(`x' + exp(`u')) * `a' * (`t_a' ^ (`a' - 1)) * exp(-exp(`x' + exp(`u')) * (`t_a' ^ `a')))) ///
                     - log(             (exp(-exp(`x'           ) * (`t_0' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_0' ^ `a')))) ///
 if `t_0' > 0  & `t_a' == `t_b'

 // [E] stock sampled, right censored
 qui replace    `lj' = log(             (exp(-exp(`x'           ) * (`t_a' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_a' ^ `a')))) ///
                     - log(             (exp(-exp(`x'           ) * (`t_0' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_0' ^ `a')))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & missing(`t_b')
 
 // [F] stock sampled, interval censored
 qui replace    `lj' = log(             (exp(-exp(`x'           ) * (`t_a' ^ `a')) - exp(-exp(`x'           ) * (`t_b' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_a' ^ `a')) - exp(-exp(`x' + exp(`u')) * (`t_b' ^ `a')))) ///
                     - log(             (exp(-exp(`x'           ) * (`t_0' ^ `a')))  ///
                           + exp(`o') * (exp(-exp(`x' + exp(`u')) * (`t_0' ^ `a')))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & !missing(`t_b') 

 mlsum `lnf' = `lj'
 
 if (`todo' == 0 | `lnf' >= .) exit // d0 ends here
 
 ********************************************
 * first derivative of log-likelihood wrt X *
 ********************************************
 
 tempvar g1j d1 

 // [A] flow sampled, uncensored
 qui gen double `g1j' = -((`t_a'^`a'*exp(`x')-1)*exp(`t_a'^`a'*exp(`x'+exp(`u')))+(`t_a'^`a'*exp(`x'+2*exp(`u')+`o')-exp(exp(`u')+`o'))*exp(`t_a'^`a'*exp(`x')))/(exp(`t_a'^`a'*exp(`x'+exp(`u')))+exp(`t_a'^`a'*exp(`x')+exp(`u')+`o')) ///
 if `t_0' == 0 & `t_a' == `t_b'

 // [B] flow sampled, right censored
 qui replace    `g1j' = -(`t_a'^`a'*exp(`x')*(exp(`t_a'^`a'*exp(`x'+exp(`u')))+exp(`t_a'^`a'*exp(`x')+exp(`u')+`o')))/(exp(`t_a'^`a'*exp(`x'+exp(`u')))+exp(`t_a'^`a'*exp(`x')+`o')) ///
 if `t_0' == 0 & `t_a' <  `t_b' & missing(`t_b')

 // [C] flow sampled, interval censored
 qui replace    `g1j' = (exp(`o')*(`t_b'^`a'*exp(-`t_b'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u'))-`t_a'^`a'*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')))+`t_b'^`a'*exp(`x'-`t_b'^`a'*exp(`x'))-`t_a'^`a'*exp(`x'-`t_a'^`a'*exp(`x')))/(exp(`o')*(exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x'))) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' != 0 & !missing(`t_b')

 // [C2] flow sampled, left censored
 qui replace    `g1j' = (`t_b'^`a'*exp(-`t_b'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`t_b'^`a'*exp(`x'-`t_b'^`a'*exp(`x')))/(exp(`o')*(1-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))-exp(-`t_b'^`a'*exp(`x'))+1) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' == 0 & !missing(`t_b')

 // [D] stock sampled, uncensored
 qui replace    `g1j' = (`a'*`t_a'^(`a'-1)*(1-`t_a'^`a'*exp(`x'+exp(`u')))*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`a'*`t_a'^(`a'-1)*(1-`t_a'^`a'*exp(`x'))*exp(`x'-`t_a'^`a'*exp(`x')))/(`a'*`t_a'^(`a'-1)*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`a'*`t_a'^(`a'-1)*exp(`x'-`t_a'^`a'*exp(`x')))-(-`t_0'^`a'*exp(-`t_0'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_0'^`a'*exp(`x'-`t_0'^`a'*exp(`x')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' == `t_b'
 
 // [E] stock sampled, right censored
 qui replace    `g1j' = (-`t_a'^`a'*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_a'^`a'*exp(`x'-`t_a'^`a'*exp(`x')))/(exp(`o'-`t_a'^`a'*exp(`x'+exp(`u')))+exp(-`t_a'^`a'*exp(`x')))-(-`t_0'^`a'*exp(-`t_0'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_0'^`a'*exp(`x'-`t_0'^`a'*exp(`x')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & missing(`t_b')
 
 // [F] stock sampled, interval censored
 qui replace    `g1j' = (exp(`o')*(`t_b'^`a'*exp(-`t_b'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u'))-`t_a'^`a'*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')))+`t_b'^`a'*exp(`x'-`t_b'^`a'*exp(`x'))-`t_a'^`a'*exp(`x'-`t_a'^`a'*exp(`x')))/(exp(`o')*(exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x')))-(-`t_0'^`a'*exp(-`t_0'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_0'^`a'*exp(`x'-`t_0'^`a'*exp(`x')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & !missing(`t_b')
 
 // Putting it all together 
 mlvecsum `lnf' `d1' = `g1j', eq(1)

 ********************************************
 * first derivative of log-likelihood wrt A *
 ********************************************

 tempvar g2j d2

 // [A] flow sampled, uncensored
 qui gen double `g2j' = -(((`t_a'^`a'*log(`t_a')*exp(`x')-log(`t_a'))*`a'-1)*exp(`t_a'^`a'*exp(`x'+exp(`u')))+((`t_a'^`a'*log(`t_a')*exp(`x'+2*exp(`u')+`o')-log(`t_a')*exp(exp(`u')+`o'))*`a'-exp(exp(`u')+`o'))*exp(`t_a'^`a'*exp(`x')))/(`a'*(exp(`t_a'^`a'*exp(`x'+exp(`u')))+exp(`t_a'^`a'*exp(`x')+exp(`u')+`o'))) ///
 if `t_0' == 0 & `t_a' == `t_b'
 
 // [B] flow sampled, right censored
 qui replace    `g2j' = -(`t_a'^`a'*log(`t_a')*exp(`x')*(exp(`t_a'^`a'*exp(`x'+exp(`u')))+exp(`t_a'^`a'*exp(`x')+exp(`u')+`o')))/(exp(`t_a'^`a'*exp(`x'+exp(`u')))+exp(`t_a'^`a'*exp(`x')+`o')) ///
 if `t_0' == 0 & `t_a' <  `t_b' & missing(`t_b')
 
 // [C] flow sampled, interval censored
 qui replace    `g2j' = (exp(`o')*(`t_b'^`a'*log(`t_b')*exp(-`t_b'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u'))-`t_a'^`a'*log(`t_a')*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')))+`t_b'^`a'*log(`t_b')*exp(`x'-`t_b'^`a'*exp(`x'))-`t_a'^`a'*log(`t_a')*exp(`x'-`t_a'^`a'*exp(`x')))/(exp(`o')*(exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x'))) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' != 0 & !missing(`t_b')
 
 // [C2] flow sampled, left censored
 qui replace    `g2j' = (`t_b'^`a'*log(`t_b')*exp(-`t_b'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`t_b'^`a'*log(`t_b')*exp(`x'-`t_b'^`a'*exp(`x')))/(exp(`o')*(1-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))-exp(-`t_b'^`a'*exp(`x'))+1) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' == 0 & !missing(`t_b')
 
 // [D] stock sampled, uncensored
 qui replace    `g2j' = (-`t_a'^(2*`a'-1)*log(`t_a')*`a'*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+2*`x'+2*exp(`u')+`o')+`t_a'^(`a'-1)*log(`t_a')*`a'*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`t_a'^(`a'-1)*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_a'^(2*`a'-1)*log(`t_a')*`a'*exp(2*`x'-`t_a'^`a'*exp(`x'))+`t_a'^(`a'-1)*log(`t_a')*`a'*exp(`x'-`t_a'^`a'*exp(`x'))+`t_a'^(`a'-1)*exp(`x'-`t_a'^`a'*exp(`x')))/(`t_a'^(`a'-1)*`a'*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`t_a'^(`a'-1)*`a'*exp(`x'-`t_a'^`a'*exp(`x')))+(`t_0'^`a'*log(`t_0')*exp(-`t_0'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')+`t_0'^`a'*log(`t_0')*exp(`x'-`t_0'^`a'*exp(`x')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' == `t_b'
 
 // [E] stock sampled, right censored
 qui replace    `g2j' = (-`t_a'^`a'*log(`t_a')*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_a'^`a'*log(`t_a')*exp(`x'-`t_a'^`a'*exp(`x')))/(exp(`o'-`t_a'^`a'*exp(`x'+exp(`u')))+exp(-`t_a'^`a'*exp(`x')))-(-`t_0'^`a'*log(`t_0')*exp(-`t_0'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_0'^`a'*log(`t_0')*exp(`x'-`t_0'^`a'*exp(`x')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & missing(`t_b')
 
 // [F] stock sampled, interval censored
 qui replace    `g2j' = (exp(`o')*(`t_b'^`a'*log(`t_b')*exp(-`t_b'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u'))-`t_a'^`a'*log(`t_a')*exp(-`t_a'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')))+`t_b'^`a'*log(`t_b')*exp(`x'-`t_b'^`a'*exp(`x'))-`t_a'^`a'*log(`t_a')*exp(`x'-`t_a'^`a'*exp(`x')))/(exp(`o')*(exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x')))-(-`t_0'^`a'*log(`t_0')*exp(-`t_0'^`a'*exp(`x'+exp(`u'))+`x'+exp(`u')+`o')-`t_0'^`a'*log(`t_0')*exp(`x'-`t_0'^`a'*exp(`x')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & !missing(`t_b')
 
 // Putting it all together 
 mlvecsum `lnf' `d2' = `g2j', eq(2)

 ********************************************
 * first derivative of log-likelihood wrt U *
 ********************************************

 tempvar g3j d3

 // [A] flow sampled, uncensored
 qui gen double `g3j' =  -((`t_a'^`a'*exp(exp(`u')+`x')-1)*exp(exp(`u')+`u'+`t_a'^`a'*exp(`x')+`o'))/(exp(`t_a'^`a'*exp(exp(`u')+`x'))+exp(exp(`u')+`t_a'^`a'*exp(`x')+`o')) ///
 if `t_0' == 0 & `t_a' == `t_b'

 // [B] flow sampled, right censored
 qui replace    `g3j' =  -(`t_a'^`a'*exp(exp(`u')+`u'+`t_a'^`a'*exp(`x')+`x'+`o'))/(exp(`t_a'^`a'*exp(exp(`u')+`x'))+exp(`t_a'^`a'*exp(`x')+`o')) ///
 if `t_0' == 0 & `t_a' <  `t_b' & missing(`t_b')
 
 // [C] flow sampled, interval censored
 qui replace    `g3j' = (exp(`o')*(`t_b'^`a'*exp(-`t_b'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x')-`t_a'^`a'*exp(-`t_a'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x')))/(exp(`o')*(exp(-`t_a'^`a'*exp(exp(`u')+`x'))-exp(-`t_b'^`a'*exp(exp(`u')+`x')))-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x'))) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' != 0 & !missing(`t_b')
 
 // [C2] flow sampled, left censored
 qui replace    `g3j' = (`t_b'^`a'*exp(-`t_b'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x'+`o'))/(exp(`o')*(1-exp(-`t_b'^`a'*exp(exp(`u')+`x')))-exp(-`t_b'^`a'*exp(`x'))+1) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' == 0 & !missing(`t_b')
 
 // [D] stock sampled, uncensored
 qui replace    `g3j' = (`a'*`t_a'^(`a'-1)*(exp(`u')-`t_a'^`a'*exp(exp(`u')+`u'+`x'))*exp(-`t_a'^`a'*exp(exp(`u')+`x')+exp(`u')+`x'+`o'))/(`a'*`t_a'^(`a'-1)*exp(-`t_a'^`a'*exp(exp(`u')+`x')+exp(`u')+`x'+`o')+`a'*`t_a'^(`a'-1)*exp(`x'-`t_a'^`a'*exp(`x')))+(`t_0'^`a'*exp(-`t_0'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x'+`o'))/(exp(`o'-`t_0'^`a'*exp(exp(`u')+`x'))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' == `t_b'
 
 // [E] stock sampled, right censored
 qui replace    `g3j' = (`t_0'^`a'*exp(-`t_0'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x'+`o'))/(exp(`o'-`t_0'^`a'*exp(exp(`u')+`x'))+exp(-`t_0'^`a'*exp(`x')))-(`t_a'^`a'*exp(-`t_a'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x'+`o'))/(exp(`o'-`t_a'^`a'*exp(exp(`u')+`x'))+exp(-`t_a'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & missing(`t_b')
 
 // [F] stock sampled, interval censored
 qui replace    `g3j' = (exp(`o')*(`t_b'^`a'*exp(-`t_b'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x')-`t_a'^`a'*exp(-`t_a'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x')))/(exp(`o')*(exp(-`t_a'^`a'*exp(exp(`u')+`x'))-exp(-`t_b'^`a'*exp(exp(`u')+`x')))-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x')))+(`t_0'^`a'*exp(-`t_0'^`a'*exp(exp(`u')+`x')+exp(`u')+`u'+`x'+`o'))/(exp(`o'-`t_0'^`a'*exp(exp(`u')+`x'))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & !missing(`t_b')
 
 // Putting it all together 
 mlvecsum `lnf' `d3' = `g3j', eq(3)

 ********************************************
 * first derivative of log-likelihood wrt O *
 ********************************************

 tempvar g4j d4
 
 // [A] flow sampled, uncensored
 qui gen double `g4j' = -((exp(`t_a'^`a'*exp(`x'+exp(`u')))-exp(`t_a'^`a'*exp(`x')+exp(`u')))*exp(`o'))/((exp(`o')+1)*(exp(`o'+`t_a'^`a'*exp(`x')+exp(`u'))+exp(`t_a'^`a'*exp(`x'+exp(`u'))))) ///
 if `t_0' == 0 & `t_a' == `t_b'

 // [B] flow sampled, right censored
 qui replace    `g4j' = exp(`o'-`t_a'^`a'*exp(`x'+exp(`u')))/(exp(`o'-`t_a'^`a'*exp(`x'+exp(`u')))+exp(-`t_a'^`a'*exp(`x')))-exp(`o')/(exp(`o')+1) ///
 if `t_0' == 0 & `t_a' <  `t_b' & missing(`t_b')
 
 // [C] flow sampled, interval censored
 qui replace    `g4j' = ((exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))*exp(`o'))/((exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))*exp(`o')-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x')))-exp(`o')/(exp(`o')+1) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' != 0 & !missing(`t_b')
 
 // [C2] flow sampled, left censored
 qui replace    `g4j' =  ((1-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))*exp(`o'))/((1-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))*exp(`o')-exp(-`t_b'^`a'*exp(`x'))+1)-exp(`o')/(exp(`o')+1) ///
 if `t_0' == 0 & `t_a' <  `t_b' & `t_a' == 0 & !missing(`t_b')

 // [D] stock sampled, uncensored
 qui replace    `g4j' = -((exp(`t_a'^`a'*exp(`x'+exp(`u'))+`t_0'^`a'*exp(`x'))-exp(`t_0'^`a'*exp(`x'+exp(`u'))+`t_a'^`a'*exp(`x')+exp(`u')))*exp(`o'))/((exp(`o'+`t_0'^`a'*exp(`x'))+exp(`t_0'^`a'*exp(`x'+exp(`u'))))*(exp(`o'+`t_a'^`a'*exp(`x')+exp(`u'))+exp(`t_a'^`a'*exp(`x'+exp(`u'))))) ///
 if `t_0' > 0  & `t_a' == `t_b'
 
 // [E] stock sampled, right censored
 qui replace    `g4j' = exp(`o'-`t_a'^`a'*exp(`x'+exp(`u')))/(exp(`o'-`t_a'^`a'*exp(`x'+exp(`u')))+exp(-`t_a'^`a'*exp(`x')))-exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & missing(`t_b')
 
 // [F] stock sampled, interval censored
 qui replace    `g4j' = ((exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))*exp(`o'))/((exp(-`t_a'^`a'*exp(`x'+exp(`u')))-exp(-`t_b'^`a'*exp(`x'+exp(`u'))))*exp(`o')-exp(-`t_b'^`a'*exp(`x'))+exp(-`t_a'^`a'*exp(`x')))-exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))/(exp(`o'-`t_0'^`a'*exp(`x'+exp(`u')))+exp(-`t_0'^`a'*exp(`x'))) ///
 if `t_0' > 0  & `t_a' <  `t_b' & !missing(`t_b')
 
 // Putting it all together 
 mlvecsum `lnf' `d4' = `g4j', eq(4)

 matrix `g' = (`d1', `d2', `d3', `d4')
end
